{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Romberg 算法计算积分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def Romberg_quadrature_algorithm(f, a, b, epsabs=1e-6):\n",
    "    \"\"\"\n",
    "    Romberg 算法计算数值积分\n",
    "    :param f: 待积分的单变量函数\n",
    "    :param a, n: 积分上下界\n",
    "    :param epsabs: 精度 default = 1e-6\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "    def cal_r():\n",
    "        def cal_c():\n",
    "            def cal_s():\n",
    "                def cal_t():\n",
    "                    n = 1\n",
    "                    t_n = (b - a) * (f(b) + f(a)) / 2\n",
    "                    while True:\n",
    "                        t_2n = 0\n",
    "                        for i in range(1, 1 + n):\n",
    "                            t_2n += f(a + (b - a) * (2 * i - 1) / (n * 2))\n",
    "                        t_2n *= (b - a) / (n * 2)\n",
    "                        t_2n += t_n * 0.5\n",
    "                        yield t_2n\n",
    "                        n += 1\n",
    "                        t_n = t_2n\n",
    "                t = cal_t()\n",
    "                while True:\n",
    "                    t1, t2 = next(t), next(t)\n",
    "                    yield 4 * t2 / 3 - t1 / 3\n",
    "            s = cal_s()\n",
    "            while True:\n",
    "                s1, s2 = next(s), next(s)\n",
    "                yield 16 * s2 / 15 - s1 / 15\n",
    "        c = cal_c()\n",
    "        while True:\n",
    "            c1, c2 = next(c), next(c)\n",
    "            yield 64 * c2 / 63 - c1 / 63\n",
    "    r = cal_r()\n",
    "    while True:\n",
    "        r1, r2 = next(r), next(r)\n",
    "        if np.abs(r2 - r1) <= epsabs:\n",
    "            return r2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 任务 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "精度为 0.1 时，圆周率的数值解为 3.141945045003288\n",
      "精度为 0.01 时，圆周率的数值解为 3.141945045003288\n",
      "精度为 0.001 时，圆周率的数值解为 3.1416767843681015\n",
      "精度为 0.0001 时，圆周率的数值解为 3.1416767843681015\n",
      "精度为 1e-05 时，圆周率的数值解为 3.1416133175296457\n",
      "精度为 1e-06 时，圆周率的数值解为 3.141597778393853\n",
      "精度为 1e-07 时，圆周率的数值解为 3.1415939298964064\n",
      "精度为 1e-08 时，圆周率的数值解为 3.141592953042991\n",
      "精度为 1e-09 时，圆周率的数值解为 3.141592720078548\n"
     ]
    }
   ],
   "source": [
    "eps_list = [10 ** i for i in range(-1, -10, -1)]\n",
    "for eps in eps_list:\n",
    "    print(\n",
    "        f\"精度为 {eps} 时，圆周率的数值解为 {Romberg_quadrature_algorithm(lambda x:4/(1+x ** 2), 0., 1., eps)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 任务2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "轨道长度：48707.439 km\n"
     ]
    }
   ],
   "source": [
    "global a, e\n",
    "h1, h2, R = 439, 2384, 6371\n",
    "c = (h2 - h1) / 2\n",
    "a = c + R + h1\n",
    "e = c / a\n",
    "\n",
    "\n",
    "def r(theta):\n",
    "    return a * (1 - e ** 2) / (1 - e * np.cos(theta))\n",
    "\n",
    "\n",
    "def f(theta):\n",
    "    return r(theta) * np.sqrt(1 + e ** 2 * np.sin(theta) ** 2 / ((1 - e * np.cos(theta)) ** 2))\n",
    "\n",
    "\n",
    "print(f\"轨道长度：{round(Romberg_quadrature_algorithm(f, -np.pi, np.pi), 3)} km\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
